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ABSTRACT 

A  procedure  for  predicting  underground  cavity  motions  resulting  from 
surface  air-blast  loading  is  developed  on  the  basis  of  a  dynamic  -  reciprocity 
principle  for  a  nonhomogeneous,  anisotropic,  linear,  viscoelastic  half 
space.  A  computer  program  for  evaluating  the  integral  equations  for  this 
medium  is  presented.  Experimental  measurements  from  the  RAINIER 
Event  of  Operation  PLUMBBOB  are  used  to  postulate  cavity  motions  arising 
from  a  surface  burst  at  the  RAINIER  site.  The  effects  of  variations  in  the 
cavity  pressure  signature,  number  of  measured  ground  motions,  and 
instrumentation  errors  are  examined.  Recommendations  are  made  for 
experimental  programs  to  obtain  needed  additional  data,  and,  in  general, 
to  ascertain  the  applicability  of  the  reciprocity  to  soils. 
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PART  I  --  SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS 


I.  INTRODUCTION 


The  central  aim  of  this  investigation  was  to  develop  a  procedure  in 
which  ground  motion  data  obtained  during  an  underground  explosion  could 
be  used  to  predict  the  underground  motion  that  would  result  if  there  was  a 
surface  explosion  at  that  site.  Part  I  summarizes  the  work  performed 
during  the  study,  while  Part  II  presents  the  details. 

During  this  study,  a  digital- computer  program  was  developed  which 
makes  it  possible  to  compute  cavity  deformations  caused  by  a  specified 
pressure  loading  on  the  surface  of  the  ground.  By  using  this  program 
the  practicality  of  the  proposed  procedure  is  demonstrated  and  no  insur¬ 
mountable  problems  are  anticipated  in  its  use. 

The  ground  motion  and  surface  pressure  pulse  data  needed  to  make 
the  cavity  motion  calculations  have  been  taken  during  numerous  tests  and 
there  are  no  difficulties  involved  in  using  these  data  in  the  digital  compu¬ 
tations.  Nevertheless,  there  still  remains  one  important  question  to  be 
answered;  that  is,  the  degree  of  applicability  of  this  procedure  to  soil 
media.  This  can  only  be  determined  when  the  necessary  data  become 
available.  Now  it  is  true  that  soil  motions  due  to  a  transient  loading,  say 
at  point  C  ,  have  been  measured,  say  at  point  S.  However,  to  test  the 
dynamic  -  reciprocity  theorem  in  a  soil,  a  reciprocal  measurement  must 
be  made.  That  is,  the  motion  at  C  due  to  a  transient  loading  at  S  should 
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be  measured.  Such  data  are  not  available  at  present,  and  suggestions  for 
a  test  program  to  gather  the  necessary  data  are  outlined  at  the  end  of 
this  report. 

The  specific  objectives  of  the  present  studies  were  formulated 
and  the  study  initiated.  These  four  objectives  were: 

« 

1)  Transform  the  general  convolution  integrals  in 
the  reciprocity  theorem  to  forms  that  are 
applicable  to  the  specified  problem  and  that 
are  suitable  for  digital  computation;  develop 

a  practical  procedure  for  data  reduction. 

2)  Investigate  supplementary  techniques  for 
relating  average  radial  displacements  to 
probable  displacement  variations  around 
the  cavity. 

3)  Evaluate  the  limitations  and  possibilities  of 
the  basic  concept  of  dynamic  reciprocity  to 

the  specified  problem;  delineate  possible  a 

extensions  to  nonlinear  and  dissipative 

continuous  systems,  and  to  other  problems  , 

in  ground  shock  propagation. 

4)  Develop  criteria  for  the  types  and  accuracy 
of  instrumentation  needed  for  field  measure¬ 
ment;  establish  reasonable  estimates  of 

the  number  and  distribution  of  measurement 
stations;  investigate  interpolation  procedures 
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for  estimating  ground  displacement 
between  stations. 

Briefly,  the  work  undertaken  to  meet  the  specific  objectives 
was  as  follows:  First,  the  dynamic- reciprocity  relation  which  is 
pertinent  to  the  problem  was  transformed  to  a  form  suitable  for  digital 
computations.  This  relation  was  programmed  and  checked  out  by  com¬ 
paring  the  computed  results  against  some  known  solutions  for  beams. 
Calculations  of  cavity  displacements  caused  by  a  postulated  nuclear 
surface  detonation  were  then  performed.  Ground-motion  velocity 
data  taken  during  the  underground  PLUMBBOB  series  of  tests  were 
used  in  these  calculations.  The  underground  cavity  pressure  which 
caused  these  ground  motions  was  not  available;  hence  a  set  of  hypo¬ 
thetical  cavity-pressure  pulse  curves  were  used.  Various  computer 
runs  were  made  to  investigate  the  effect  on  the  computed  cavity 
deformations  caused  by  various  cavity  pressure  pulse  curves.  The 
influence  of  errors  in  the  ground- station  velocity  measurements  was 
also  investigated. 

Since  the  cavity  displacement  calculations  provide  only  the  average 
radial  displacement,  the  relation  between  the  average,  distortional,  and 
rigid-body  radial  displacements  was  investigated  next.  Published  theo¬ 
retical  results  from  a  related  problem  were  used  to  estimate  the  relative 
magnitude  of  these  radial  displacements.  The  related  problem  was  solved 
by  Baron  and  Parnes  (Reference  6)  and  dealt  with  the  response  of  a  cylin¬ 
drical  cavity  during  the  passage  of  a  plane  compressional  wave.  They 
determined  the  radial  and  tangential  velocities  of  the  cavity  boundary  as  a 
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function  of  time  for  several  points  around  the  periphery  of  the  cavity  in 
nondimens ional  variables.  From  these  published  results  the  average, 
rigid  body,  and  distortiona.l  displacement  time-histories  for  several  points 
around  the  cavity  were  determined.  This  problem  differed  from  the  spheri¬ 
cal  cavity  problem  considered  in  this  study  in  two  respects.  First,  the 
medium  was  assumed  to  be  homogeneous,  isotropic,  nondissipative,  and 
linearly  elastic;  and  second,  the  cavity  was  cylindrical  rather  than  spheri¬ 
cal.  In  an  effort  to  remove  the  objection  of  a  different  geometry,  the 
analogous  problem  for  a  spherical  cavity  was  examined  and  its  solution 
formulated.  However,  the  calculations  of  the  cavity  displacements  re¬ 
quired  the  development  of  another  computer  program.  Since  no  numerical 
computations  were  performed  the  derivation  of  this  formal  solution  to  the 
spherical  cavity  problem  was  not  included  in  this  report  but  has  been  docu¬ 
mented  in  STL  Report  EM  13-7. 

II.  DIGITAL  PROGRAM  DEVELOPMENT 

A  special  form  of  the  general  dynamic-reciprocity  theorem,  (derived 
in  the  Appendix),  which  has  been  transformed  and  programmed  on  the  IBM 
7  090  so  that  calculations  of  cavity  motions  could  be  made  is 

» 

[  f  pc(t  "  T)uc(T)dScdT  =  f  f  Pg(t  ~  t)u  <T)dS  dt  (1) 

Jo  Js  Jo  J  S  * 

c  s 

where  p  and  u  represent  pressure  and  displacement  time-histories  and  sub¬ 
scripts  c  and  s  refer  to  cavity  and  surface.  This  convolution  integral  re¬ 
lation  is  integrated  over  the  cavity  area  S^  and  the  surface  area  Sg  as  well 
as  over  time.  In  the  applications  of  this  integral  relation  it  is  assumed 
that  the  surface  responses  ug  caused  by  the  cavity  pressure  pc  have  been 
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measured  and  the  cavity  deformation  caused  by  the  specified  surface 
pressure  pg  is  to  be  determined.  Upon  examination  of  the  type  of  ground 
motion  data  taken  during  underground  explosions,  it  appeared  that  ground 
velocity  data  ii  or  acceleration  data  u  were  more  readily  obtainable  than 
ground  displacement  data.  It  was  decided,  therefore,  to  transform  ( 1 )  into 
a  form  where  the  digital  computations  of  the  unknown  cavity  deformation 
uc(t)  depended  upon  the  ground  velocities  ug,  and  the  two  pressure  pulses 
p^  and  pg.  In  the  computational  scheme  adopted,  the  space  integration 
over  the  ground  surface  was  replaced  by  a  summation  over  n  stations  by 
assuming  that  the  variations  of  the  integrands  over  each  station  were  solely 
time  variations.  Clearly,  this  requires  that  the  spacing  between  stations 
be  small  enough  so  that  the  variations  of  pg  and  ug  from  station  to  station 
are  not  too  large.  From  a  computational  standpoint  this  does  not  place  any 
serious  restrictions  upon  the  cavity  displacement  calculations  since  the 
computer  is  capable  of  handling  a  large  number  of  stations.  Generally, 
the  upper  limit  to  the  number  of  surface  stations  will  be  governed  by  the 
number  of  velocity  time-histories  that  can  be  recorded  during  the  under¬ 
ground  detonation.  However,  from  the  data  to  be  presented  presently,  it 
appears  that  the  number  of  surface  stations  required  are  within  the  capa¬ 
bilities  of  the  usual  instrumentation  procedures. 

Before  the  digital  program  could  be  used  to  make  cavity  deformation 
calculations  it  was  necessary  to  obtain  surface  and  cavity  pressures  and 
ground  velocity  time-histories.  In  the  calculations  it  was  assumed  that 
the  surface  pressures  were  caused  by  a  surface  nuclear  explosion.  Once 
a  given  yield  and  ground  zero  were  selected,  Brode's  theoretical  pressure 
pulse  was  determined  for  each  of  the  n  stations  (References  2  and  14). 
Actually,  this  was  found  to  be  a  rather  tedious  task  when  a  large  number 
of  stations  used,  and  therefore,  an  alternate  method  for  introducing 
the  pressure  pulse  data  was  devised.  The  digital  program  was  modified 
so  that  an  arbitrary  pressure  pulse  curve  at  each  station  can  be  input. 

Or  alternatively,  the  weapon  yield  and  the  distance  from  ground  zero  to 
each  of  the  n  stations  are  specified,  then  the  computer  generates  the 
appropriate  Erode  pressure  pulse  curve. 
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III.  CAVITY  DISPLACEMENT  CALCULATIONS 


A  search  for  data  on  ground  motion  caused  by  underground  detonation 
showed  that  a  few  measurements  had  been  made  on  a  large  number  of  tests, 
but  there  were  no  tests  in  which  surface  motions  had  been  measured  at  a 
large  number  of  stations  along  several  lines  radiating  out  from  surface 
ground  zero.  The  most  complete  ground  motion  data  found  were  ground- 
velocity  data  taken  during  the  Rainier  test  in  the  PLUMBBOB  series. 

During  this  test,  ground  velocities  were  measured  at  seven  different  sta¬ 
tions  ranging  from  ground  zero  out  to  a  range  of  2340  feet.  The  locations, 
of  the  underground  detonation  and  respective  ground  stations  for  these  data 
are  shown  in  figure  1.  The  stations  were  very  nearly  on  a  single  radial  line 
emanating  from  surface  ground  zero. 


STATION  1  2  3  4  5  6  7 


Figure  1.  Location  of  Surface  Stations  with  Respect 
to- Unde  rground  Detonation — Rainier  Test 

Since  no  ground- motion  data  were  available  to  indicate  the  nonisotropic 
nature  of  the  medium  where  these  data  were  taken,  it  was  tacitly  assumed 
that  the  medium  was  isotropic,  and  therefore,  these  measured  motions 
were  identical  with  those  along  any  other  radial  line. 

The  ground- motion  velocity  time  -  historie  s  at  each  of  these  stations  are 
shown  in  figures  7  and  8  in  Part  II.  Unfortunately,  no  data  were  available 
on  the  cavity  pressure  p^  that  produced  the  ground  velocities,  and  there¬ 
fore,  a  theoretical  pressure  pulse  was  used.  For  the  calculations,  it  was 
assumed  that  the  cavity  pressure  could  be  expressed  in  the  form 

p  (t)  -  Ae'at  +  Be'Pt  (2) 

c 
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The  peak  pressure,  i.e.  ,  A  +  B,  was  selected  from  the  theoretical  peak 
pressure  curve  for  various  cavity  radii  presented  in  figure  8  of  Reference 
14  for  an  assumed  effective  radius  of  about  10  feet.  This  corresponds  to 
a  pressure  of  about  100,000  atmospheres.  The  remaining  three  constants 
in  equation  (2)  were  arbitrarily  selected.  These  parameters  were  changed 
in  some  of  the  computations  to  investigate  the  influence  of  variations  in  the 
peak  cavity  pressure  and  decay  rate. 

Numerous  computations  were  made  using  the  ground  velocity  data  and 
cavity  pressure  pulses  for  a  surface  nuclear  explosion  having  its  ground 
zero  at  one  of  two  locations.  The  first  location  was  taken  as  directly  above 
the  underground  cavity,  and  therefore,  coincident  with  ground  zero  for  the 
underground  explosion.  This  case  is  referred  to  as  the  direct-hit  example. 
The  second  location  for  ground  zero  was  1880  feet  away  from  the  ground 
zero  for  the  underground  explosion.  This  case  is  referred  to  as  the  non- 
direct-hit  example.  The  conclusions  drawn  from  these  computations,  and 
to  be  presented  subsequently,  are  not  meant  to  be  interpreted  as  conclusions 
that  are  valid  in  general,  for  clearly,  they  are  strongly  dependent  upon  the 
data  used.  Nevertheless,  the  calculations  are  probably  indicative  of  a 
typical  situation  and  should  provide  a  guide  to  the  influence  of  the  various 
parameters.  In  the  calculations  it  was  assumed  that  the  medium  was  iso¬ 
tropic  and,  therefore,  the  measured  velocity  time-history  at  any  station  is 
the  same  for  any  other  station  at  the  same  distance  from  surface  ground 
zero. 

A.  INFLUENCE  OF  CAVITY  PRESSURE  PULSE  SIGNATURE 

The  firshparameter  studied  is  the  influence  of  various  cavity  pressure 
pulse  time-histories  upon  the  computed  cavity  displacements.  This  is  an 
important  parameter  to  investigate  for  there  is  difficulty  in  measuring  this 
pressure  pulse  and  theoretical  computations  are  strongly  dependent  upon 
the  properties  of  the  explosive  material,  cavity  size,  and  the  simplifying 
assumptions.  For  example,  see  the  theoretical  pressure  pulse  curves 
computed  in  Reference  3. 

Figure  3  shows  the  computed  cavity  displacements  for  the  direct-hit 
example  for  the  three  cavity  pressure  pulses  having  peak  pressures  of 
100,  000,  81,  625,  and  61, 000  atmospheres  shown  in  figure  2.  The  pressure 
pulse  curves  were  selected  by  initially  taking  the  pressure  pulse  to  be 
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TIME  (SECONDS) 

Figure  Z.  Assumed  Cavity  Pressures  During  the  Rainier  Shot 
in  the  PLUMBBOB  Series 
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CAVITY  AVERAGE  RADIAL  DISPLACEMENT,  u  (FEET) 


0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6 

TIME  (SECONDS) 


Figure  3.  Cavity  Displacements  for  Direct-Hit  Example 
Based  on  Cavity  Pressures  of  Figure  2. 
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-0. 25t 


(3) 


p  (t)  =  75,  000  e'50t  +  25,  000  e 
c 

and  then  selecting  new  values  of  A,  B,  a,  and  p  subject  to  the  condition  that 

pc(t>  at 

is  a  constant.  Note  that  in  figure  3,  the  peak  cavity  displacements  are  1.2, 
0.75,  0.  55  ft  for  the  three  pressure  pulses,  p^,  p  ^ ,  Pcg>  having  peak 

pressures  of  100,000,  81,625,  and  6  1 ,  000  atmospheres  respectively .  These 
data  show  that  if  there  were  indeed  as  much  uncertainty  in  the  pressure 
pulse  curve  as  between  these  extremes,  there  would  be  an  uncertainty  of 
approximately  2  to  1  in  the  computed  maximum  cavity  displacement.  This 
suggests  that  when  this  procedure  is  applied  to  predict  the  cavity  displace¬ 
ment  caused  by  a  surface  pressure,  efforts  should  be  taken  to  properly 
define  the  cavity  pressures. 

From  figure  2  note  that  p  ^  is  a  larger  pressure  over  the  first  1.  6 
seconds  except  for  the  initial  instant,  while  figure  3  shows  that  this  cavity 
pressure  pulse  curve  leads  to  the  smallest  predicted  cavity  displacement. 
This  indicates  that  if  one  uses  too  large  a  cavity  pressure,  the  computed 
cavity  displacements  will  be  too  small.  Obviously,  it  is  not  necessary  to 
make  cavity  deformation  calculations  to  note  this  fact  for  it  is  evident 
from  equation  (1).  If  one  assumes  pg  and  ug  are  known  functions  of  time, 
then  the  right  side  of  equation  (1)  is  some  fixed  time  function.  Thus,  if  on 
the  left  side  p^  is  too  large,  then  u^  will  be  too  small  accordingly.  Physi¬ 
cally,  what  this  means  is  that  if  a  small  cavity  pressure  produces  a  large 
surface  motion,  then  similarly  a  small  surface  pressure  produces  a  large 
cavity  motion.  These  comments  suggest  that  if  there  is  any  uncertainty  in 
the  cavity  pressure  pulse,  it  is  conservative  to  use  a  cavity  pressure  curve 
that  is  too  small. 

B.  INFLUENCE  OF  NUMBER  OF  SURFACE  MEASUREMENTS 

Several  calculations  were  made  to  investigate  the  effect  of  the 
number  of  surface  stations  used  in  the  computations  upon  the  computed 
cavity  deformations.  One  such  calculation  was  made  by  using  every  other 
surface  station  for  the  same  direct-hit  example  used  above.  Adopting  the 
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station  numbering  system  shown  in  figure  1  for  the  direct-hit  example,  calcu¬ 
lations  were  made  by  using  the  four  stations  1,  3,  5,  and  7  and  the 
results  were  compared  with  those  when  all  seven  stations  were  used.  The 
results  of  this  computation  are  shown  in  figure  4.  This  figure  shows  that 
the  predicted  cavity  deformations  were  larger  when  four  stations  were  used. 
This  is  primarily  a  consequence  of  the  manner  of  adjusting  the  surface 
areas  in  the  two  sets  of  computations.  In  the  four  station  computations,  the 
area  associated  with  station  1  was  larger  than  when  all  seven  stations  were 
used.  Now,  since  the  surface  pressure  and  velocity  were  the  largest  at 
station  1,  the  contribution  of  station  1  to  the  total  cavity  deformation  is 
large.  Therefore,  increasing  the  area  of  station  1  increases  the  predicted 
cavity  response.  This  example  provides  an  illustration  of  a  case  where  the 
spatial  variation  of  the  time  functions  was  too  large  between  stations,  at 
least,  for  those  stations  closest  to  ground  zero.  This  does  not  necessarily 
infer  that  four  stations  are  insufficient,  but  it  does  show  that  if  fewer  sta¬ 
tions  are  used  and  the  maximum  pressure  and  velocity  occurring  within  the 
area  associated  with  the  respective  station  are  assumed  to  act  over  the  en¬ 
tire  area,  then  a  conservative  estimate  of  cavity  motion  will  be  obtained. 

In  an  effort  to  provide  additional  information  regarding  the  number  of 
surface  stations  and  the  effect  of  instrumentation  errors  at  these  stations, 
some  influence  functions  were  computed  for  both  the  direct-hit  and  nondirect- 
hit  examples.  These  influence  functions  were  defined  as  the  cavity  response 
caused  by  the  surface  pressures  acting  on  each  station  individually.  In  the 
results  summarized  here,  the  cavity  pressure  pulse  given  in  equation  (3) 
was  used. 

The  station  numbering  system  shown  in  figure  1  is  again  used  for  the 
direct-hit  example.  Figure  5  shows  the  maximum  cavity  displacement 
produced  by  the  pressures  acting  on  each  of  the  seven  stations.  Probably 
it  should  be  mentioned  that  these  maximum  cavity  displacements  are  pro¬ 
portional  to  the  area  lumped  at  each  surface  station  as  well  as  to  the  meas¬ 
ured  ground  velocity  and  prescribed  overpressure.  For  example,  the  area 
lumped  at  station  7  was  about  three  times  that  lumped  at  station  6.  This 
larger  area  associated  with  station  7  caused  the  influence  of  station  7  to 
be  greater  than  station  6.  However,  on  an  influence  per  unit  area  basis, 
the  influence  of  successive  stations  would  clearly  decrease  as  the  distance 
from  the  cavity  to  the  station  increased. 
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Figure  4.  Comparison  of  Cavity  Displacements 
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Figure  5.  Influence  Coefficients  —  Maximum  Cavity  Displacements 
Caused  by  Pressures  Acting  on  Individual  Stations 
for  the  Direct-Hit  Example 
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Figure  5  is  developed  from  the  cavity  motion  data  presented  in 
figures  14  and  15  of  Part  II.  Note  that  figure  5  shows  that- the  absolute 
value  of  the  maximum  cavity  displacement  is  1.21  feet  while  the  maxi¬ 
mum  value  from  the  first  station  alone  is  0.78  foot.  Since  the  maximum 
cavity  response  caused  by  each  surface  station  does  not  occur  at  the  same 
instant  of  time,  the  maximum  cavity  response  is  not  obtained  by  adding  the 
maximum  values  shown  in  figure  5.  From  figure  3,  it  is  noted  that  the 
maximum  cavity  displacement  occurs  at  0.76  second  for  the  pressure 
pulse  curve  p  j.  Figure  6  shows  the  individual  contribution  of  each 
station  to  this  maximum  cavity  displacement  at  this  instant  of  time.  In 
this  figure,  it  is  noted  that  the  contribution  from  stations  4,  5,  6,  and  7 
are  small  compared  to  the  first  three  stations.  In  fact,  if  only  the  first 
three  stations  are  used,  the  computed  cavity  displacement  is  -1.  14  feet, 
rather  than  the  value  of  -1.21  feet  when  all  seven  stations  are  used.  In 
such  a  case,  the  computed  error  using  the  first  three  stations  would  be 
about  6  percent.  Recall  that  when  the  four  stations  1,  3,  5,  and  7  were 
used,  the  error  was  approximately  70  percent.  Thus,  reasonable  pre¬ 
dictions  of  cavity  displacements  can  be  made  with  a  limited  number  of 
stations;  however,  the  stations  must  be  selected  judiciously  as  indicated 
later  on. 

For  the  nondirect-hit  example,  the  station  numbering  system  is  shown 
in  figure  7.  More  stations  are  required  for  the  nondirect-hit  example  than  for 
the  direct-hit  example  because  of  the  loss  of  some  symmetry.  Figure  8 
shows  the  maximum  cavity  displacement  caused  by  each  individual  station  when 
considered  separately.  Indicated  in  this  figure  is  a  resultant  maximum 
cavity  displacement  of  0.26  foot  when  the  properly  phased  contribution 
from  all  stations  is  summed.  Figure  9  shows  the  actual  contribution  of 
each  station  to  this  maximum  cavity  displacement.  Since  the  cavity  motion 
caused  by  pressures  acting  on  some  stations  begins  after  the  maximum 
cavity  displacement  is  reached,  their  contributions  to  the  maximum  are 
zero.  From  figure  9,  it  is  noted  that  if  only  the  six  stations,  1,  la,  lb, 

2,  2a,  and  2b  had  been  used  in  the  analysis,  the  predicted  maximum  cavity 
displacement  would  have  been  -0.268  foot,  which  is  very  close  to  the 
value  of  -0.260  foot,  when  all  surface  stations  are  considered. 
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In  both  the  direct-hit  and  nondirect-hit  calculations,  a  fairly  accurate 
cavity  displacement  can  be  predicted  using  a  limited  number  of  surface 
stations,  provided  that  they  are  properly  selected.  The  most  important 
stations  are  those  stations  where  the  major  portion  of  the  energy  from  the 
nuclear  surface  detonation  enters  the  ground,  i.e.,  in  the  vicinity  of  ground 
zero.  In  the  calculations  shown  here  for  an  isotropic  medium  and  for  the 
direct-hit  example  reasonably  accurate  cavity  displacements  were  predicted 
using  only  three  stations  when  full  advantage  of  the  problem  symmetry  is 
used.  For  the  nondirect-hit  example,  some  symmetry  was  lost  and  fairly 
accurate  cavity  displacements  were  predicted  using  six  stations. 

C.  INFLUENCE  OF  ERRORS  IN  SURFACE  MEASUREMENTS 

Errors  in  measured  surface  velocities  cause  a  corresponding  error  in 
the  computed  cavity  deformations.  For  example,  since  a  linear  system  is 
assumed,  say  a  10  percent  error  in  every  measured  ground  velocity  pro¬ 
duces  a  10  percent  error  in  the  computed  cavity  deformation.  However, 
instrumentation  errors  may  be  randomly  distributed  so  that  the  cancelling 
effects  prod  me  u  smaller  error.  Figure  ^  can  be  used  to  estimate  the 
error  in  the  predicted  cavity  deformation  Cor  a  specified  error  in  either  the 
surface  pressure  or  ground  velocity  at  any  station.  Clearly,  the  largest 
error  in  predicted  cavity  deformation  results  from  an  error  in  the  meas¬ 
urements  taken  at  station  1.  Figure  5  shows  that  the  ratio  of  the  contri¬ 
bution  from  station  1  to  the  contribution  from  all  the  stations  is  approxi¬ 
mately  0.67/1. 21  =0.  bs.  Therefore,  if  say  die  surface  velocity  data  for 
station  1  had  a  constant' 10  percent  error,  this  would  produce  an  error  of 
a.  5  percent  in  the  predicted  cavity  deformation.  Computations  of  this  type 
led  to  the  conclusion  that  errors  in  the  computed  cavity  deformation  caused 
by  instrumental  errors  in  the  data  taken  at  the  ground  surface  stations  are 
secondary  when  compared  with  the  errors  caused  by  the  theoretical  approxi¬ 
mation  to  the  actual  pressure  pulse  curve,  and  the  errors  resulting  by  re¬ 
placing  the  space-time  integration  of  the  convolution  integral  over  the  sur¬ 
face  by  summation  of  time  integration  at  each  of  the  n  surface  stations. 
These  conclusions  indicate  that  the  typical  accuracy  obtainable  in  ground 
motion  and  pressure  measurements  is  entirely  satisfactory.  The  accuracy 
of  cavity  motion  predictions  can  be  improved  more  effectively  by  increas¬ 
ing  the  number  of  surface  measurements  rather  than  by  improving  the 
accuracy  of  the  measurements  made  at  a  given  set  stations. 


D.  GROUND  MOTIONS  BETWEEN  MEASURED  POINTS 


The  possibility  was  examined  of  interpolating  between  stations  to 
determine  the  ground  motions  at  some  intermediate  stations.  Here  the 
hope  was  that  fewer  ground  motion  measurements  could  be  made  and  the 
accuracy  of  the  cavity  computations  improved  by  increasing  the  number  of 
stations  used  in  the  analysis.  However,  it  was  soon  apparent  that  no  ac¬ 
curacy  could  be  gained  by  this  method  and  one  might  just  as  well  use  those 
stations  for  which  data  were  available. 


IV.  CAVITY  DISPEACEMENTS  — AVERAGE,  DISTORTION, 
AND  TRANSLATION 


Assume  ground  surface  velocity  data  are  obtained  from  a  symmetric 
pressure  pulse  having  a  zero  resultant  force  acting  on  an  underground 
cavity.  If  these  data  are  used  to  compute  cavity  deformations  caused  by 
a  pressure  loading  on  the  ground  surface,  only  average  cavity  deformations 
can  be  determined  by  the  dynamic  reciprocity  procedure.  The  fact  that  the 
procedure  provides  no  information  regarding  the  cavity  distortion  and  rigid- 
body  displacement  is  not  a  limitation  of  the  dynamic- reciprocity  procedure, 
but  rather  emphasizes  the  fact  that  the  cavity  pressure  pulse,  which  pro¬ 
vided  the  free-surface  velocity  data,  was  symmetric  and  had  a  zero  re¬ 
sultant  force.  In  fact,  if  a  nonsymmetric  pressure  pulse  having  a  resultant 
force  can  be  produced  inside  the  cavity,  and  if  the  spatial  and  time  distri¬ 
bution  of  this  pressure  pulse  is  known,  then  these  pressure  pulse  data  and 
the  associated  ground  velocity  data  can  be  used  to  determine  the  total  cavity 
displacement  that  would  result  from  a  specified  ground- surface  pressure 
loading.  It  is  quite  likely,  however,  that  when  this  dynamic- reciprocity 
procedure  is  applied  the  cavity  pressure  pulse  will  be  symmetric  and  have 
a  zero  resultant  force .  For  this  reason,  theoretical  computations  were 
conducted  to  see  how  the  average  displacement  was  related  to  the  distortion 
and  rigid-body  displacement.  In  this  study  use  was  made  of  some  results 
presented  by  Baron  and  Parnes  in  Reference  6. 


Baron  and  Parnes  determine  the  motion  of  the  cylindrical  cavity  in  the 
infinite  elastic  medium  illustrated  in  figure  10  during  the  passage  of  a  plane 
compressional  wave.  Presented  in  their  work  are  plots  of  the  nondimen- 
sional  radial  velocities  [w/(trc  /pfj  and  tangential  velocity  jjv / ( cr c  /p)J  ,  as 
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the  cavity,  where  \v  is  the  radial  velocity,  v  the  tangential  velocity,  cr  the 
stress  intensity  of  the  wave,  p  the  shear  modulus  of  the  medium,  the 
propagation  velocity  of  the  dilatational  waves,  a  the  radius  of  the  cylindri¬ 
cal  cavity  'and  t  the  time. 


Figure  10.  Plane  Compre s s ional  Wave  Engulfing  a 
Cylindrical  Cavity  in  an  Infinite  Medium 

The  detailed  work  presented  in  Part  II  consisted  of  integrating  the 
cavity  velocity  components  to  obtain  the  displacement  components  (0,t) 
and  ( 0 ,  t) .  Next,  the  displacements  were  expanded  into  a  Fourier  series 

to  determine  individual  contributions  of  the  average,  distortional,  and  rigid- 
body  motion  to  the  total  radial  displacement.  This  was  done  by  representing 
the  cavity  displacements  by  the  equations 

co 

U  (0,  t)  =  a  (t)  +  y  a  (t)  cos  n0 
c  o  n 

n=  1 

co 

V  (0,  t)  =  y  b  (t)  sin  n0 
c  n 

n=  1 

so  that  the  component  displacements  are  given  by 
average  radial  displacement  =  aQ(t) 

a  j.  (t)  -  b  1  (t) 

rigid-body  radial  displacement  =  - 2 -  cos  ® 

aj(t)  -  bj (t) 

distortional  radial  displacement  =  U  (9,  t)  -  aQ(t)  -  - j -  cos  ® 
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The  results  of  these  computations  showed  that  the  maximum  displace¬ 
ments  occurred  at  0  =  0.  Figure  11  shows  a  plot  of  the  three  displacements 
at  0  =  0,  as  a  function  of.  time.  From  this  figure  it  is  noted  that  the  average 
displacement  and  distortional  displacement  are  roughly  of  the  same  magni¬ 
tude.  This  suggests  that,  if  only  the  average  radial  displacement  is  com¬ 
puted,  the  total  elastic  displacement  will  be  obtained  approximately  by 
multiplying  by  a  factor  of  two.  Note  that  after  the  wave  front  has  propagated 
one  diameter  beyond  the  cavity  (nondimensional  time  equal  to  two),  the 
cavity  translational  motion  increases  nearly  at  a  constant  rate.  This  con¬ 
tinuing  translational  motion  is  a  consequence  of  the  assumption  made  in  this 
problem  that  the  compressive  stress  is  continually  applied.  Thus,  as  the 
compressive  wave  front  advances  with  a  constant  velocity,  the  entire  medium 
moves  with  a  constant  velocity  without  acceleration.  However,  such  rigid- 
body  translations  without  distortion  of  the  cavity  shape  do  not  cause  any 
stress  variations. 


V.  SUMMARY 

During  this  study,  a  dynamic- reciprocity  procedure  was  devised  and  a 
digital-computer  program  developed  whereby  ground-motion  data  taken 
during  an  explosion  in  an  underground  cavity  could  be  used  to  predict  cavity 
motions  which  would  be  caused  by  a  specified  pressure  pulse  loading  on  the 
ground  surface.  This  digital  program  was  so  arranged  that  the  data  re¬ 
quired  to  compute  cavity  motion  are  the  pressure  pulse  time-history  inside 
the  cavity,  the  corresponding  normal  velocities  of  the  ground  surface,  and 
the  prescribed  pressure  loading  on  the  ground  surface.  Given  these  data, 
the  computer  evaluates  the  average  radial  cavity  deformation  time-history 
that  the  specified  surface  pressure  loading  induces.  If  a  nuclear  detonation 
is  prescribed  on  the  surface,  the  yield  of  the  weapon  and  the  distance  from 
ground  zero  to  each  of  several  surface  stations  can  be  specified  and  the 
computer  will  generate  the  approximate  Brode  pressure  pulse  curve  as  is 
necessary  during  the  cavity  calculations.  This  procedure  is  rigorously 
valid  in  nonisotropic,  linear  viscoelastic  media,  and,  therefore,  can  be 
used  for  linear  dissipative  media.  Even  though  the  program  was  written 
to  use  ground  velocities,  acceleration  measurements  can  be  used. 
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The  required  number  of  ground  motion  measurements  that  must  be 
recorded  so  as  to  apply  the  dynamic-reciprocity  procedure  is  very  much 
dependent  upon  the  nature  of  the  medium.  For  example,  if  the  medium  is 
isotropic,  measurements  along  one  radial  line  are  sufficient.  Whereas,  if 
the  medium  is  nonisotropic  measurements  should  be  made  along  three  or 
four  radial  lines.  The  number  of  measurements  required  along  any  radial 
line  should  be  at  least  3,  but  from  4  to  7  or  more  are  recommended.  With 
the  data  used  in  this  study,  it  was  noted  that  the  error  in  the  calculated 
average  cavity  displacement  was  about  SO  percent  using  only  one  surface 
measurement  and  about  13  percent  using  three  surface  measurements, 
when  compared  to  the  calculated  average  cavity  displacement  using  seven 
surface  measurements. 

Cavity  distortions  can  be  estimated  from  the  average  cavity  displace¬ 
ment  using  the  results  from  the  plane  wave  diffraction  problem.  Roughly, 
the  maximum  distortion  is  about  equal  to  the  average  displacement.  An 
alternative  and  more  satisfactory  method  of  estimating  the  cavity  distor¬ 
tion  and  also  translation  can  be  made  using  experimentally  determined 
transfer  fum-ions. 


VI.  RECOMMENDATIONS 

This  study  demonstrated  that  a  practical  procedure  utilizing  dynamic- 
reciprocity  relations  can  be  used  to  predict  the  average  cavity  displace¬ 
ments  caused  by  a  prescribed  pressure  heading  acting  on  the  surface,  pro¬ 
viding  experimental  measurements  of  ground  motions  caused  by  a  known 
average  pressure  acting  inside  the  cavity  are  available.  Although  the 
average  radial  displacement  is  vital  to  the  designer  of  underground  struc¬ 
tures,  the  distortional  and  translational  displacements  are  also  important. 
Theoretically,  dynamic  reciprocity  provides  a  means  for  predicting  total 
cavity  motion,  including  translation  and  distortion  as  well  as  the  average 
displacement  caused  by  a  hypothetical  pressure  acting  on  the  surface  of  the 
ground.  However,  such  predictions  can  be  made  only  if  one  knows  the  spa¬ 
tial  distribution  of  the  pressure  time-histories  acting  around  the  cavity 
while  the  ground  motion  data  are  being  collected.  If  only  the  average  pres¬ 
sure  time-history  acting  in  the  cavity  is  known,  then  only  the  average 
radial  displacement  caused  by  a  prescribed  surface  loading  can  be  predicted. 
Generally  the  pressure  distribution  around  the  cavity  will  be  nearly  uniform 
unless  some  clever  experimental  procedure  is  adopted. 


2  4 


It  is  clear  that  requiring  detailed  cavity  pressure  loading  imposes  a 
severe  constraint  upon  the  test  procedure  followed  when  collecting  the  ground 
motion  data.  For  this  reason,  it  is  desirable  to  remove  this  requirement 
by  adopting  an  alternate  procedure.  A  more  practical  procedure  is  to  use 
the  combined  experimental-theoretical  procedure  which  only  requires  the 
average  cavity  pressure,  for  predicting  the  average  radial  displacement 
caused  by  a  prescribed  ground  pressure  loading.  Then,  using  the  average 
displacement,  along  with  some  additional  experimental  data,  estimate  the 
distortion  and  translation  of  the  cavity.  To  use  this  alternate  procedure, 
data  are  required  which  show  the  relation  among  the  cavity  average,  dis- 
tortional,  and  rigid-body  displacements.  Such  data  can  be  obtained  most 
conveniently  on  a  linear  medium  using  scaled  models. 

In  the  appendix,  the  dynamic  reciprocity  theorem  is  derived  and  shown 
to  be  valid  for  a  nonhomogeneous  viscoelastic  medium  provided  the  stress- 
strain  relationship  is  linear.  At  the  present  time  proof  does  not  exist  of 
the  validity  of  the  reciprocity  relation  in  a  nonlinear  medium.  Neither  is 
there  a  proof  to  the  contrary.  Lamb,  in  Reference  16,  states  on  the  basis 
of  Lagrange's  results  that  the  restriction  to  linearity- may  not  be  necessary, 
but  he  did  not  elaborate  on  this  point,  and  a  proof  is  lacking. 

Clearly,  in  many  applications  of  dynamic  reciprocity  one  must  deal 
with  nonlinear  media.  For  example,  while  rocks  may  behave  elastically, 
soils  generally  do  not.  It  is  worthwhile,  therefore,  to  find  out  experi¬ 
mentally  the  degree  of  applicability  of  the  uynamic- reciprocity  relation¬ 
ship  to  such  media. 

It  is  recommended  that  a  simple  and  definitive  experimental  program 
be  initiated  that  will  fulfill  the  following  two  objectives: 

1)  Gather  sufficient  test  data  to  determine  the  relation  among 
cavity  translation,  distortion,  and  average  displacement 
in  a  linear  medium. 

2)  Gather  the  necessary  data  to  determine  the  degree  of  appli¬ 
cability  of  dynamic  reciprocity  in  a  nonlinear  medium  such 
as  a  soil. 

The  first  objective  will  provide  the  transfer-function  type  of  data  required 
for  estimating  cavity  distortions  and  rigid-body  translations  from  the  aver¬ 
age  displacement  for  a  linear  medium.  The  second  objective  will  make  it 
possible  to  estimate  the  errors  in  the  predicted  cavity  displacements  caused 
by  nonlinear  effects. 
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PART  II  -  ANALYTICAL  RESULTS 


I.  TRANSFORMATION  OF  THE  RECIPROCITY  RELATION 

A  special  form  of  the  dynamic  -  reciprocity  relation,  equation  (16)  in 
the  Appendix,  which  is  suitable  for  our  purposes  is  given  by 


f  rt  r 

rt 

Is  Jo  Pc(Xlt-T)uc(x’T)  dTdsc  =JS 
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(1 

where  p  is  the  pressure,  u  is  the  displacement,  S  is  the  area,  x  repre¬ 
sents  the  space  coordinates,  and  t  and  t  are  time  coordinates.  The 
subscripts  s  and  c  represent  surface  and  cavity  respectively.  The 
problem  becomes  one  of  solving  equation  (1)  for  the  unknown  .cavity  dis¬ 
placement  u  caused  by  a  specified  surface  pressure  loading  p  when  the 
normal  surface  motion  ug  has  been  measured  due  to  a  known  cavity  pres¬ 
sure  p^.  To  solve  this  integral  equation,  several  simplifying 
assumptions  are  introduced.  First,  it  is  realized  that  the  ground  surface 
measurements  will  only  be  known  at  a  finite  number  of  surface  stations 
so  that  space  integration  over  the  ground  surface  must  be  replaced  by  a 
summation  over  the  n  measured  stations.  It  is  assumed  that  at  each  of 
these  n  stations,  both  the  pressure  and  ground  motion  are  the  same  over 
the  entire  station  area  at  any  instant  of  time.  It  is  further  assumed  that 
the  internal  pressure  p^  acting  on  the  cavity  wall  is  uniformally  distri¬ 
buted  and  simply  varies  with  time.  Under  this  assumption  u^  becomes 
the  average  radial  cavity  displacement  produced  by  the  ground  surface 
pressure  loading  p  .  It  is  also  assumed  that  the  variation  of  the  cavity 
area  with  time  can  be  neglected  and  that  this  area  can  be  taken  as  a  con¬ 
stant.  These  assumptions  permit  equation  (1)  to  be  rewritten  as 
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where  Sc  represents  the  cavity  area,  Sgn  represents  the  ground  surface 
area  of  the  ntd  station  and  the  summation  is  to  extend  over  the  n  ground 
surface  stations.  To  cast  equation  (2)  into  a  form  suitable  for 
programming  on  the  IBM  7090,  it  is  necessary  that  several  algebraic 
manipulations  be  performed. 
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Due  to  the  symmetry  property  of  the  convolution  integral,  equation  (2) 
may  be  rewritten  as 
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After  equation  (3)  is  differentiated  with  respect  to  time  t,  it  becomes 
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where  dots  denote  differentiation  with  respect  to  the  independent  variable 

of  the  functions.  It  is  assumed  that  the  initial  displacements  at  the  ground 

surface  u  (0)  are  zero  so  that  the  last  term  in  equation  (4)  vanishes, 
sn 

Then,  applying  the  symmetry  property  of  the  convolution  integral  and 
solving  equation  (4)  for  the  average  radial  displacement  of  the  cavity 
boundary  gives 
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In  equations  ( 4)  and  ( 5)  is  the  velocity  of  the  ground  surface  due  to  an 

underground  explosion  and  p<_  is  the  rate  of  the  pressure  change  at  the 
cavity  boundary  during  the  explosion. 

By  examining  equation  (5)  it  is  observed  that  the  first  integral  contains 

the  known  functions  p  and  u  whereas  the  second  integral  contains  the 
rsn  sn 

known  function  p^  and  the  unknown  function  u^.  It  is  convenient  to  separate 
each  of  these  two  integrals  into  two  parts.  The  first  integral  will  be 
separated  into  two  parts  to  save  computational  time  and  the  second  inte¬ 
gral  into  two  parts  to  express  entirely  in  terms  of  known  functions, 
thereby  avoiding  the  iteration  procedure  that  would  be  required  to  solve 
for  u^  from  equation  (5). 
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In  the  first  integral  in  equation  (5),  it  is  recognized  that  the  over¬ 
pressure  p  will  be  zero  until  the  pressure  arrives  at  station  n  at  the 
tifne  t^  .  Therefore,  the  lower  limit  in  this  integral  can  be  changed 
from  t  =  0  to  t  -  t^.  Now,  in  general,  the  pressure  p  decays  rapidly 
after  the  arrival  of  the  peak  overpressure  and  then  is  followed  by  a  more 
slowly  decaying  pulse,  especially  when  these  pressures  are  caused  by  a 
nuclear  detonation.  Computational  time  can  be  reduced  by  avoiding  the 
small  integration  step  required  to  evaluate  the  contribution  from  the  fast 
decay  portion  by  analytically  integrating  the  fast  decay  portion  and  then 
numerically  integrating  the  remaining  portion.  By  letting  e  represent 
the  duration  of  the  fast  decay  part  of  p  ,  the  integral  is  replaced  by  an 

integral  covering  the  period  from  t  =  t.  to  t  =  t,  +  e  and  another  inte- 
'  In  In  n 

gral  for  the  period  from  t  =  t,  +  e  to  t  =  t. 

In  n 

The  second  integral  is  separated  into  two  parts,  covering  the  period 
from  t  =  0  to  t  =  At,  and  from  t  =  At  to  t  =  t  where  At  represents  the 
time  interval  selected  for  evaluating  the  cavity  displacements.  In  this 
manner,  the-  integral  from  At  to  t  contains  only  known  functions,  since 
when  u^  is  being  determined  at  time  t  it  will  be  known  for  all  times  up 
to  (t  -  At).  However,  the  integral  from  0  to  At  contains  the  time  portion 
of  u^  that  is  being  determined,  so  that  this  integral  will  require  further 
examination.  After  these  manipulations,  equation  (5)  is  rewritten 
as 
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For  sufficiently  small  €  ,  the  value  of  u  can  be  considered  constant 
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during  the  period  t^  to  (t^  +  e  )  and  the  approximation 
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is  made  during  this  period.  For  this  period,  the  pressure  p  is 
approximated  by  a  double  exponential  decay  equation 
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where  A  ,  B  ,  n  ,  B  are  constants.  If  equations  (7)  and  (8)  are 
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substituted  into  equation  (6),  the  first  integral  in  equation  (b)  becomes 
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which  upon  integration  gives 
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In  a  similar  manner,  it  is  assumed  that  the  cavity  pressure  p^  can 
be  represented  by  a  double  exponential  decay  equation  during  the  time 
interval  from  0  to  At.  Thus,  p^  is  taken  as 
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p  (  t)  =  A  e  C+Be  C  .  0  £  t  S  At  (11) 
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where  attain  A  .  B  .  a  .  and  B  are  constants.  Previously,  it  was  pointed 
=  c  c  c  r  c 

out  that  when  u  is  being  determined  at  time  t.  u  will  be  known  for  the 
entire  period  from  0  to  (t  -  At)  so  that  the  last  integral  in  equation  (6) 
involves  an  integration  of  known  functions.  However,  the  third  integral 
in  equation  (6)  concerns  that  period  of  time  for  which  u^  is  not  yet  known. 
To  circumvent  this  difficulty  and  express  this  integral  in  terms  of  the 
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values  of  at  the  beginning  and  end  of  time  interval  0  -  t  -  At,  another 
approximation  is  introduced.  The  approximation  is  made  that  the  varia¬ 
tion  in  u  (t)  which  represents  the  function  u^(t  -  t)  in  the  interval  0  -•  t 
At  can  be  expressed  by  the  parabola 


u  (t)  =  aT^  +  br  +  c  , 
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subjected  to  the  conditions: 
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where  u  ^  =  the  radial  displacement  of  the  cavity  at  the  end  of  the  previous 
time  interval  (t  -  At)  and  u^.  =  the  radial  velocity  of  the  cavity  at  the  end 
of  the  previous  time  interval  (t  -  At).  By  imposing  these  conditions, 
equation  (13)  becomes 
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and  the  function  u  (t  -  r)  is 
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Then  the  third  integral  of  equation  (6)  may  be  written  as 
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where  the  procedure  of  integration  by  parts  has  been  used. 

Finally,  by  introducing  equations  (10),  (11),  and(18)  into  equation  ( 6) , 
the  cavity  displacement  u^  becomes 
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II.  DIGITAL  PROGRAM 


An  IBM  7090  digital  program  was  developed  to  solve  equation  (19) 
for  the  cavity  radial  motion  u^.  To  use  this  program,  the  velocity  and 
pressure  data  must  first  be  introduced  into  the  computer.  Some  of  the 
methods  of  introducing  these  data  and  several  other  salient  points  will 
be  discussed  next  before  presenting  computed  results. 

A.  GROUND  VELOCITIES  u 

sn 

The  digital  program  is  so  arranged  that  ground  velocity  data  meas¬ 
ured  normal  to  the  ground  surface  and  caused  by  an  underground  explosion 
are  used.  These  data  are  input  into  the  computer  in  tabular  form,  that  is, 
the  measured  velocities  at  a  discrete  set  of  times.  If  ground  acceleration 
or  displacement  data  are  recorded  during  the  underground  explosion,  it 
will  be  necessary  to  integrate  or  differentiate  such  data  to  derive  the 
corresponding  velocity  data  before  they  are  introduced  into  the  computer. 
Clearly,  acceleration  data  are  preferable  to  displacement  data  to  avoid 
any  loss  of  accuracy  through  the  differentiation  process. 

B.  SURFACE  PRESSURES  p 

rsn 

The  surface  pressures  acting  on  each  of  the  n  stations  can  be  intro¬ 
duced  into  the  computer  in  one  of  two  ways.  The  first  method  provides 
a  means  of  handling  those  cases  where  the  pressure  pulse  signature  is 
completely  arbitrary.  The  second  metnod  is  restricted  to  a  nuclear  sur¬ 
face  detonation.  This  second  method  of  introducing  the  pressure  pulse 
data  for  a  nuclear  detonation  was  developed  to  simplify  the  input  procedure. 
This  seemed  to  be  a  worthwhile  simplification  because  it  was  anticipated 
that  most  problems  to  be  considered  will  be  those  in  which  the  cavity 
motion  must  be  determined  for  a  nuclear  explosion  on  the  surface. 

1  .  Arbitrary  Pressure  Pulse 

The  procedure  for  introducing  the  pressure  pulse  acting  at 
each  of  the  n  stations  is  identical  to  that  used  for  the  ground  velocities. 
That  is,  a  tabular  listing  of  the  discrete  set  of  times  and  pressures  is 
stored  in  the  computer.  A  slight  variation  of  this  procedure  can  be  fol¬ 
lowed  when  the  pressure  pulse  decays  rapidly  following  the  peak  pressure. 
In  this  variation,  the  fast-decay  portion  is  approximated  by  the  double 
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decaying-exponential  pulse  as  expressed  by  equation  (8),  and  a  tabular 
listing  lor  the  remaining  time  interval.  The  analytical  pre s sure  pulse 
is  generated  internally  in  the  computer  as  required  during  the  computa¬ 
tions,  once  the  constants  A  ,  B  ,  n  ,  (3  ,  and  t,  are  prescribed  for 

sn  sn  sn  '  sn  in  r 

each  station. 


2.  Nuclear  Pressure  Pulse 

The  overpressure  acting  at  each  of  the  n  stations  caused 
by  a  nuclear  surface  detonation  is  also  generated  within  the  computer 
once  the  weapon  yield  and  the  distance  from  ground  zero  to  each  of  the 
n  stations  are  specified.  This  generated  overpressure  pulse  is  an  approx¬ 
imation  to  the  theoretical  overpressures  derived  by  Erode  (Reference  2). 
These  overpressure  pulses  were  approximated  by  the  double  decaying- 
exponential  equation 
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for  the  subsequent  time  period,  where  «  is  .-et  to  0.  Oi  second.  The 

values  of  A  ,  B  .a  ,6  ,  t,  ,  and  the  corresponding  primed  param- 
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eters  are  internally  evaluated  in  the  computer  by  fitting  the  above  two 
equations  to  Erode' s  data  for  the  specified  range  from  ground  zero  and  for 
a  nominal  yield. 

Once  the  weapon  yield  W  and  the  distances  of  each 
station  n  from  ground  zero  are  specified,  the  computer  evaluates  the 
pressures  p  and  arrival  times  t^  by  the  following  steps.  First  the 
computer  evaluates  the  equivalent  distance  R^,  for  the  nominal  yield 
weapon  W^_  at  which  the  same  overpressures  would  be  observed  using 
the  equation  (Reference  1  1). 
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The  value  of  is  then  used  to  evaluate  the  constants  A  ,  B  ,  a  ,6  , 

E  sn  sn  sn  rsn 

the  corresponding  primed  constants  and  the  arrival  times  Uj^)n  for  the 
nominal  yield  weapon.  Finally,  the  overpressures  Psn>  acting  on  each 
station,  are  computed  according  to  equations  (20)  and  (21)  and  the  arrival 
times  at  these  stations  from  the  equation 
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C.  CAVITY  PRESSURE  p 


The  cavity  pressure  time-history  can  be  introduced  into  the  computer 
in  one  of  two  ways.  Either: 


1) 


Prescribe  the  constants  Ac', 


Br 


and  the  computer 


evaluates  the  pressure  at  any  time  t  from  the  equation 
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or 


2)  Prescribe  the  constants  Ac,  Bc,  ac,  |3C,  and  the  computer 

evaluates  the  pressure  for  the  first  time  interval  0  to  At  using 
equation  (11).  The  cavity  pressures  following  this  period 
are  prescribed  by  a  tabular  listing  of  discrete  pressure  and 
time  values. 

D.  INTEGRATION  STEP  SIZE 

The  integration  step  size  can  be  specified  in  one  of  two  ways  in  the 
Adams  -  Moulton,  Runge-Kutt.a  method  used  for  numerically  integrating 
the  integrals  in  equation  (19).  Either  a  constant  step  size  can  be  pre¬ 
scribed,  or  limits  placed  on  the  truncation  error  allowed  in  any  integra¬ 
tion  step.  In  the  latter  method,  a  variable  step  size  is  used  with  the 
computer  selecting  its  own  step  size  by  comparing  the  error  introduced 
at  each  step  with  the  prescribed  allowable  error,  and  then,  adjusting  the 
step  size  to  keep  the  error  within  the  specified  limits. 
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III.  PROGRAM  CHECKOUT 


A.  BEAM  EXAMPLE 

The  programming  of  equation  (1),  and  also  the  steps  leading  to 
equation  (19)  from  the  reciprocity  relation  (1)  were  verified  by  computa¬ 
tions  for  problems  in  which  exact  solutions  were  known  or  obtainable. 
One  such  problem  used  for  this  purpose  was  the  response  at  some  point 
of  a  simply  supported  beam  due  to  transient  loadings  applied  at  three 
other  points. 


>1  s2  s3 


rthn  B 


Figure  1.  Beam  Geometry 

Figure  1  shows  the  beam  used  for  this  check  where  the  three  loading 
points  are  si,  s 2,  s3  and  the  response  is  determined  at  point  c.  In  order 
that  this  problem  can  bear  some  resemblance  to  the  half  space  problem 
to  be  considered  later,  the  loadings  at  the  three  points  approximate  a 
pressure  pulse  traveling  across  the  beam  with  different  arrival  times 
and  reduced  peak  pressures  at  successive  points  along  the  beam.  The 
loadings  selected  are  shown  in  figure  Zand  are  expressed  mathematically  as 

Pgl^*-)  =  Pll^*  "  tjj)  -  P12P^  *  ^12^  *  P13P^  '  ^13)  "  ^14^^  ” 

Ps2^  =  ^l^1  '  ^l^'  P22H^  *  ^22  ^  “  P23H^  ’  l23^  ’  P24P^  '  t24^  ,M25) 

Ps3^  =  P31H^  '  l31^  “  P32H^  *  C32 ^  '  P33H^  "  l33*  "  P34H^  '  l34t 

where  H(t)  is  the  Heaviside  unit  step  function. 
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LOADI 


It  is  found  by  analysis  that  these  loadings  produce  a  displacement  at  point  c 
given  by 


oo 


u 


I 


1 


i-  1  m  =  1  Mco 


,  mirR 
mim  n 

sm  — —  sin  — j - 


n  1 


1  -  cos  u>  (t 
m 


nl 


nZ 


1  -  cos  wm(t-tn2)]H(t-tn2)  -  Pn3 


1  -  COS  CO  (t  -  t  ^  ) 

m  n3 


1  -  cos  com(t  -  tn4)|  H(t  -  t 


I 


n4' 


H(t-tn3) 


(26) 


Now,  to  apply  equation  (19)  to  find  u  ,  the  velocity  at  the  three 
stations  must  be  determined  for  a  known  loading  at  point  c.  The  loading 
chosen  to  be  applied  at  point  c  is  shown  in  figure  3 


Figure  3.  Loading  at  Point  c 
and  is  expressed  mathematically  as 

pc(t)  =  PclH(t  -  tcl)  -Pc2Hlt  -  tc2)  -  Pc3H(t  -  tc3) 

The  velocities  at  points  s  caused  by  the  loading  (equation  (27)J  are 
found  analytically  to  be 


(27) 


mTrR 


UJ  .  V  1 1 1  M  a.  v  , 

i,„(*l  =  T-  ,  HT  sin  T11  si"  — -  pc  1 v<*  -  *c  1 : >»<•  -  ‘c  r 1 

m=  1  m  \ 

-  Pc2  Sin  “m(t  -  ‘cZ^1  -  tcZ)  "  Pc3  sin  wm(t  ‘  lc3)H(t  ’  ‘cSY  (28) 
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for  n  =  1,  2,  3.  These  velocities  were  evaluated  for  the  modes  m  =  1  to  7 
and  are  shown  in  figure  4.  The  velocity  time-history  data  and  the  pressure 
loadings  pgn  and  p^  were  introduced  into  the  computer  and  uc  determined 
with  a  constant  integration  step  size  of  0.  05  second. 

For  these  computations,  the  following  values  were  used: 


El 

1 

R2 

0.  7 

M63 

2  5it4 

St 

wn 

2 

m 

~  ~ 

R3 

i 

0.  8 

M 

=  30  lb/in 

‘cl 

0  sec 

Pcl 

=  14  lb 

h 

St 

=  0.  25 

fcc2  = 

2.  0  sec 

Pc2 

=  8  lb 

R1 

St 

=  0.  55 

lc3  = 

4.  5  sec 

Pc3 

=  6  lb 

The  results  of  these  calculations  are  shown  in  figure  5  where  the  value 

of  the  analytical  solution  of  was  determined  from  equation  (26)  by 

summing  the  terms  from  m  =  1  to  m  -  7.  The  higher  mode  contributions 

for  m  =-7  for  both  u  and  u  were  negligible  compared  to  the  first  seven 
c  sn 

mode  contributions.  Figure  5  shows  that  the  computer  solution  and  the 
analytical  solution  are  quite  close  with  the  maximum  difference  about 
5  percent.  These  results  suggested  that  the  coding  and  the  manipulation 
of  equation  (1)  into  the  form  given  by  equation  (19)  were  correct,  and 
therefore,  the  digital  program  was  then  used  to  investigate  the  cavity 
displacements  caused  by  surface  pressure  loadings. 

B.  INTEGRATION  STEP  SIZE 

Several  beam  deflection  calculations  were  made  to  investigate  the 
effect  of  the  integration  step  size  used  in  the  digital  computations.  The 
step  size  selected  in  any  given  computation  depends  on  the  data  to  be 
integrated  and  the  accuracy  demanded.  The  data  used  in  these  calculations 
are  the  loading  curves  shown  in  figure  2  and  the  velocity  curves  shown  in 
figure  4..  It  is  clear  that  when  the  loading  and  velocity  data  vary  rapidly, 

f 

the  integration  step  size  must  be  reduced  accordingly  provided  these  rapid 
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Figure  4.  Beam 
Due  to 


VERTICAL  BEAM  DEFLECTION  AT  POINT  c  (INCHES) 


TIME  (SECONDS) 


Figure  5.  Comparison  of  Analytical  Solution  and  Computer 
Solution  for  the  Beam  Problem 
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fluctuations  are  large  enough  to  significantly  affect  the  final  results.  By 
examining  figure  4,  it  is  observed  that  the  primary  velocity  curve  has  a 
rise  time  of  about  five  seconds  while  the  individual  fluctuations  have  a 
minimum  rise  time  of  about  0.  1  second.  The  discontinuities  in  the  load¬ 
ing  curve  shown  in  figure  2  also  give  rise  to  higher  frequency  fluctuations 
(this  would  be  evident  if  the  loading  curve  is  replaced  by  its  Fourier  series 
expansion).  Several  calculations  were  made  using  various  integration  step 
sizes.  Figure  6,  which  shows  th6  results  for  integration  step  sizes  of 
0.  05  and  0.  5  second,  indicates  that  the  difference  between  the  two  pre¬ 
dicted  deflections  at  t  =  5.  5  seconds  is  about  10  percent.  Figure  5  shows 
that  the  deflection  computed  by  using  0.05  second  step  size,  agrees  fairly 
well  with  the  exact  deflection,  so  that  the  integration  step  size  of  0.  05 
second  is  entirely  satisfactory,  and  even  the  0.  5  second  step  size  is 
probably  adequate  in  this  problem. 

IV.  CAVITY  DISPLACEMENT  CALCULATIONS 

To  compute  the  average  cavity  radial  displacement  caused  by  a  pre¬ 
scribed  pressure  acting  on  the  ground  surface,  it  is  necessary  to  know 
the  ground  motions  caused  by  a  known  underground  cavity  pressure  load¬ 
ing.  The  most  nearly  complete  ground  motion  data  that  could  be  found 
were  ground  velocity  measurements  obtained  during  the  Rainier  shot  in 
the  PLUMBBOB  series  (Reference  22).  These  data  consist  of  seven 
surface  measurements  that  are  nearly  on  one  radial  line  extending  from 
surface  ground  zero  out  to  a  distance  of  2340  feet.  Figures  7  and  8  show 
these  data  for  the  stations  at  distances  of  0,  450,  650,  900,  1270,  1730, 
and  2340  feet  from  surface  ground  zero.  The  cavity  pressures  that  gave 
rise  to  these  ground  motions  were  not  recorded  during  the  test  and  there¬ 
fore,  a  pressure  pulse  must  be  assumed.  Reference  14  indicates  that 
the  cavity  radius  before  the  shot  was  3  feet  and  after  the  explosion,  the 
crushed  radios  before  collapse  was  55  feet.  For  the  calculations 
which  follow  an  effective  radius  of  10  feet  will  be  assumed.  For 
this  radius,  the  cavity  area  is  1200  square  feet.  Figure  8  in 
Reference  14  shows  the  theoretical  peak  pressure  inside  the  cavity  as 
a  function  of  cavity  radius.  For  the  cavity  radius  of  10  feet,  the  peak 
pressure  is  about  100,000  atmospheres.  In  the  following  calculations, 
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VERTICAL  BEAM  DEFLECTION  AT  c  (I 


TIME  (SECONDS) 


Figure  6.  Comparison  of  Calculated  Beam  Deflection  Using 
Two  Integration  Step  Sizes 
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Figure  8.  Ground  Velocity  Data  for  the  Rainier  Shot  in  the  PLUMBBOB 
Series  for  Ranges  900,  1270,  1730  and  2340  ft 
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it  is  assumed  that  the  pressure  pulse  has  a  zero  rise  time  to  the  maximum 
pressure  and  then  decays  exponentially  as  described  by  equation  (24).  One 
pressure  pulse  used  for  many  of  the  calculations  to  follow  was  taken  as 

p  =  75,  000  e‘50t  +  25,  000  e'0,  25t  (29) 

c 


With  these  ground  velocity  data  and  cavity  pressure  pulse,  the  average 
cavity  response  was  computed  for  a  nuclear  surface  detonation  at  one  of 
two  locations.  These  two  cases  are  called  the  direct-hit  and  nondirect- 
hit  examples. 

A.  DIRECT-HIT  EXAMPLE 

For  the  direct-hit  example,  it  is  assumed  that  a  nuclear  detonation 
occurs  on  the  surface  with  its  ground  zero  coinciding  with  the  surface 
ground  zero  for  the  underground  detonation.  Brode's  theoretical  over¬ 
pressures  (Reference  2)  are  used  for  the  overpressures  acting  on  each 
of  the  ground  stations,  and  are  axisymmetr ical  about  the  ground  zero 
location. 

1 .  Ground  Surface  Stations 

It  is  assumed  that  the  medium  is  isotropic  and  homogeneous  so 
that  the  recorded  velocities  during  the  underground  explosion  would  have 
been  the  same  if  they  had  been  measured  along  any  other  radial  line  from 
surface  ground  zero.  Since  both  the  ground  pressures  p^  and  ground 
velocities  ug  are  symmetric  about  ground  zero,  the  surface  areas  used 
in  the  solution  of  equation  (19)  are  the  ring  shaped  areas  around  ground 
zero  indicated  in  figure  9.  These  surface  areas  are  listed  in  table  1. 

Table  1.  Station  Areas  in  the  Direct-Hit 
Example  Using  7  Stations 


Station 

n 


Area 

ft2 


1 

2 

3 

4 

5 

6 
7 


159.  4 
791, 300 
936, 600 
1, 81 1,  000 
3,  370, 000 
5,  941,  000 
14,  330,  000 
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2 .  Ground  Overpressures 


Since  Brode's  theoretical  overpressure  is  infinite  at  ground  zero, 
the  first  ground  station  was  arbitrarily  moved  out  to  a  point  150  feet  from 
ground  zero.  The  overpressure  used  at  this  first  station  is  that  corre¬ 
sponding  to  Brode's  overpressure  at  a  range  of  150  feet,  whereas  the 
ground  velocity  used  there  is  the  actual  velocity  recorded  at  surface 
ground  zero.  The  overpressure  data  for  the  fast  decay  portion  were  com¬ 
puted  from  equation  (20)  using  the  values  listed  in  table  2.  For  the  re¬ 
maining  period,  the  overpressures  were  input  into  computer  in  tabular 
form  by  reading  overpressure  and  time  values  from  the  overpressures 
shown  in  figure  10. 

Table  2.  Overpressure  Constants  Used  in  equation  (20) 


Station 

n 

A 

sn 

atmos 

B 

sn 

atmos 

a 

Sn-  1 
sec 

PSn_i 

sec 

lln 

sec 

i 

97,  500 

2500 

.6000 

250 

0.  01 

2 

2,  bOO 

800 

1500 

195 

0.  007 

3 

400 

200 

520 

200 

0. 0035 

4 

320 

80 

420 

25 

0. 0045 

5 

120 

10 

120 

20 

0.  01 

6 

40 

19 

70 

10 

0.  02 

7 

17 

7.  4 

31 

5 

0.  076 

3.  Influence  of  Cavity 

Pressure 

Signature 

Since  there  is  considerable  uncertainty  in  the  cavity  pressure 
pulse,  calculations  were  made  using  various  pressure  pulses  so  that  the 
influence  of  the  cavity  pressure  signature  could  be  observed.  The  basic 
pressure  pulse  is  taken  as  that  expressed  in  equation  ( <2. 9 ) .  Other  pulses 
are  selected  arbitrarily  subject  to  the  constraint  condition. 


dT  = 


constant 


Figure  12  shows  the  computed  average  cavity  displacements  caused  by  the 
nuclear  explosion  on  Station  1  for  the  five  cavity  pressures  shown  in 
figure  1 1 . 
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OVERPRESSURE  (ATMOSPHERES) 


TIME  (SECONDS) 


Figure  10.  Overpressures  Acting  on  the  Seven  Ground  Stations 
for  the  Direct-Hit  Example 


48 


TIME  (SECONDS) 

Figure  11.  Typical  Cavity  Pressures  Used  in  the  Calculations 
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4.  Influence  of  Ground  Surface  Stations 


In  order  to  estimate  the  effect  of  using  a  fewer  number  of  stations 
while  computing  cavity  displacements,  calculations  were  made  using  every 
other  station  shown  in  figure  9,  i.  e.  ,  stations  1,  3,  5,  and  7.  The  adjusted 
surface  areas  for  this  case  are  shown  in  table  3. 


Table  3.  Station  Areas  for  the  Direct-Hit 
Example  Using  Four  Stations 

Area 
103  ft2 


Station 

n 


1 

3 

5 

7 


332 
2,  560 
7,  340 
15,  700 


The  results  of  these  calculations  are  shown  in  figure  13  where  the  p^  ^ 
curve  given  by  equation  (29)  was  again  used.  In  figure  13  it  is  noted  that 
this  four-station  lumping  predicts  a  larger  cavity  motion,  This  is  pri¬ 
marily  a  consequence  of  the  increased  area  used  for  station  1  in  the  four- 
station  computations. 

Figures  14  and  15  show  the  results  of  another  set  of  calculations 
made  to  investigate  the  influence  of  each  surface  station.  These  curves 
were  computed  by  using  all  seven  stations  shown  in  figure  9.  Each  curve 
on  figures  14  and  15  shows  the  cavity  motion  caused  by  the  overpressures 
acting  on  each  of  the  seven  stations  individually.  For  example,  figure  14 
shows  that  the  maximum  cavity  displacement  caused  by  the  overpressure 
acting  on  station  1  alone  is  -0.  75  foot,  whereas  figure  12  shows  the  maxi¬ 
mum  to  be  -1.21  feet  when  the  properly  phased  contributions  from  all 
seven  stations  are  considered. 

B.  NONDIRECT-HIT  EXAMPLE 

The  nondir ect-hit  example  is  defined  as  one  in  which  the  location  of 
the  surface  nuclear  detonation  does  not  coincide  with  the  surface  ground 
zero  for  the  underground  explosion.  In  this  example,  it  is  assumed  that 
a  surface  burst  having  the  same  yield  as  in  the  direct-hit  example  takes 
place  at  a  point  1880  feet  from  surface  ground  zero.  Once  again,  Brode's 
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Figure  14.  Influence  Coefficients  —  Cavity  Displacements  Caused  by 
Surface  Pressure  Loadings  Acting  on  Stations  1,  2,  and  3 
for  the  Direct- Hit  Example 
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Influence  Coefficients  —  Cavity  Displacements  Caused  by 
Surface  Pressure  Loadings  Acting  on  Stations  4,  5,  6, 
and  7  for  the  Direct-  Hit  Example 


theoretical  overpressures  are  used  for  the  p  data  and  the  Rainier 
ground  velocities  for  the  u  data.  The  cavity  pressure  signature  used 
is  given  by  equation  (29)  and  denoted  by  p  . 

For  these  calculations,  the  38  stations  shown  in  figure  16  are  used. 
Note  that  even  though  the  surface  overpressure  has  symmetry  about  ground 
zero  for  the  surface  detonation  and  the  ground  velocities  have  symmetry 
about  the  surface  ground  zero  for  the  underground  explosion,  the  products 
of  the  overpressures  and.  velocities  only  have  symmetry  about  the  line 
passing  through  the  two  ground  zero  locations.  Thus,  when  full  advantage 
is  taken  of  this  symmetry,  the  cavity  motion  caused  by  the  surface  pres¬ 
sures  acting  on  all  38  stations  can  be  determined  by  using  the  26  stations 
1,  la,  lb,...,  20,  and  2  1 .  However,  this  is  considerably  more  than  was 
required  for  the  direct-hit  example.  Table  4  lists  the  geometrical  data 
required  in  these  calculations.  Shown  in  the  table  are  the  areas  of  each 
of  the  different  26  stations,  the  distance  from  surface  ground  zero  for 
the  underground  explosion  to  the  respective  station,  and  the  distance  from 
the  surface  burst  to  each  station. 

The  cavity-displacement  calculations  made  using  the  above  data  are 
summarized  in  figures  17  to  2  1  in  the  form  of  influence  curves.  That 
is,  each  curve  shows  the  cavity  motion  caused  by  the  overpressure  acting 
on  each  of  the  stations  individually.  The  total  cavity  motion  caused  by 
the  pressures  acting  on  all  the  stations  is  obtained  from  the  sum  of  these 
individual  curves. 

These  figures  can  be  used  to  compare  the  relative  importance  of 
those  stations  where  the  surface  pressures  are  the  largest  with  those 
where  the  ground  velocities  are  the  largest.  For  example,  compare  the 
effect  of  station  2  where  the  surface  explosion  acts  with  that  of  the  four 
stations  around  surface  ground  zero  for  the  underground  explosion.  Even 
though  table  4  shows  that  the  area  of  station  2  is  approximately  equal  to 
the  total  area  of  the  four  stations,  figures  17,  18,  and  19  show  that  the 
effect  of  station  2  is  about  twice  as  large  as  the  combined  effect  of  the  four 
stations  surrounding  surface  ground  zero.  Thus  in  this  example  the  larg¬ 
est  contribution  to  the  cavity  displacements  comes  from  the  stations 
around  ground  zero  where  the  ground  overpressures  are  the  largest. 
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Figure  16.  Lumped  Stations  for  the  Nondirect- Hit  Example 
X  =  surface  burst  ground  zero 
(J  =  underground  explosion  ground  zero 
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Table  4. 

Geometrical  Parameters  for  the  Nondirect 

•Hit  Example 

Station 

No. 

Station 

Area 
ft2  x  103 

Distance  From 
Underground 

Explosion 

Ground  Zero  (ft) 

Distance  From 
Surface  Burst 
Ground  Zero  (ft) 

1 

398 

2340 

460 

la 

597 

2340 

660 

lb 

995 

2340 

1270 

2 

165 

1730 

150 

2a 

248 

1730 

415 

2b 

413 

1730 

1030 

3 

281 

1270 

610 

3a 

281 

1270 

1010 

4 

453 

900 

980 

5 

234 

650 

1230 

6 

198 

450 

1430 

7 

39.  8 

150 

1730 

8 

39.  8 

150 

2030 

9 

198 

450 

2330 

10 

234 

650 

2530 

1  1 

453 

900 

2780 

12 

843 

1270 

3150 

13 

1485 

1730 

3610 

14 

3582 

2340 

4220 

15 

39.  8 

150 

1880 

16 

198 

450 

1930 

17 

234 

650 

1990 

18 

453 

900 

2080 

19 

843 

1270 

2270 

20 

1485 

1730 

2560 

21 

3582 

2340 

3010 
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AVERAGE  RADIAL  DISPLACEMENT,  u  (FEET) 


0.04 


) 


TIME  (SECONDS) 


Figure  17.  Influence  Coefficients — Cavity  Displacements  Caused  by 
Surface  Pressure  Loadings  Acting  on  Stations  1,  ia,  and 
lb  for  the  Nondirect-Hit  Example 
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Figure  18.  Influence  Coefficients  —  Cavity  Displacements  Caused  by 
Surface  Pressure  Loadings  Acting  on  Stations  2,  2a, 
and  2b  for  the  Nondirect- Hit  Example 
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Figure  19.  Influence  Coefficients  —  Cavity  Displacements  Caused  by 
Surface  Pressure  Loadings  Acting  on  Stations  3.  3a, 

4,  5,  6,  and  7  for  the  Nondirect-Hit  Example 
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AVERAGE  RADIAL  DISPLACEMENT,  u  ,  (FEET) 
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Figure  20. 


0.02 


-0.02 


I 

1 

- , - , - - - 1 - 

-XX _ _ _ 

1  - 

1 

n 

1 

1 

— i — i — ' — | — ; — i 

1  ! 

1 

: 

12 

_ i _ * _ _ i _ i _ 

1 

Influence  Coefficients  —  Cavity  Displacements  Caused  by 
Surface  Pressure  Loadings  Acting  on  Stations  8  to  14 
for  the  Nondirect-Hit  Example 
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V.  COMPARISON  OF  CAVITY  AVERAGE 
DISPLACEMENT  AND  DISTORTION 


Even  though  the  average  cavity  displacement  caused  by  a  specified 
nuclear  burst  on  the  surface  is  of  vital  interest,  cavity  distortions  are 
also  important.  The  numerical  procedure  used  in  the  examples  only 
predicts  the  average  cavity  displacement.  This  section  examines  a 
theoretical  problem  in  which  the  average  displacement  and  distortion 
are  computed  to  investigate  the  relative  magnitude's  of  each  of  these 
types  of  displacements.  The  calculations  of  average  displacement  and 
distortion  and  also  rigid-body  translation,  are  determined  from  some 
published  results  of  Baron  and  Parnes  (Reference  6).  Baron  and  Parnes 
find  the  motion  of  a  cylindrical  cavity  during. the  passage  of  a  plane  com- 
pressional  wave.  They  assume  that  the  medium  is  homogeneous,  iso¬ 
tropic,  and  elastic  with  a  Poisson's  ratio  of  1/4  which  corresponds  to 
a  granite-rock  medium.  The  plane  compre s sional  wave  front  has  a  step 
distribution  in  time  and  propagates  at  the  constant  velocity  of  the  dilata- 

tional  waves  c  in  a  direction  normal  to  the  generators  of  the  cylindrical 
P 

cavity.  Figure  22  shows  the  geometrical  relations  and  the  points  on  the 
cylinder  surface  where  the  radial  displacements  are  computed. 


The  average  and  distortional  radial  displacements  of  points  on  the 
cylindrical  surface  are  derived  from  the  data  shown  in  Reference  6. 
Figures  7  and  8  of  Reference  6  show  the  radial  and  tangential  velocities 
as  functions  of  time  in  terms  of  nondimensional  coordinates  (wp/a-c^  and 
vp/crc  versus  c  t/2a)  for  the  twelve  points  around  the  cavity  0^  =  0, 
e2  =  ±22.  5°,  03  =  ±45°,  04  =  ±67.  5°,  05  =  ±90°,  e6  =  ±135°,  and  0?  =  180°. 
In  these  plots  w  and  v  represent  the  radial  and  tangential  velocities  re¬ 
spectively,  a  the  stress  intensity  of  the  incoming  wave,  c^  the  propagation 
velocity  of  the  incoming  waves,  p  the  shear  modulus  of  the  medium, 
t  the  time,  and  a  the  cavity  radius.  The  radial  velocity  data  have  been 
reproduced  and  are  shown  in  figure  23  of  this  report. 

Let  Uc  and  represent  the  nondimensional  radial  and  tangential 
displacements  and  let  these  displacements  be  represented  by  the  Fourier 
expansions 


oo 

U  (0,  t)  =  a  (t)  +  Y  a  (t)  cos  n0 
c  o  ^ .  n 

n=  1 


(30) 
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Figure  ZZ.  Geometry  of  Baron  and  Parnes  Cylindrical 
Cavity  Problem 
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(31) 


oo 

V  (0,  t)  =  Y  b  (t)  sin  n0 
c  n 

n=  1 

In  this  notation,  the  nondimensional  average  and  distortional  radial  dis¬ 
placements  and  rigid-body  translations  in  the  radial  direction  can  be 
written  as 

average  radial  displacement  =  aQ(t)  (32) 

aj(t)  -  b  1  (t) 

rigid-body  radial  translation  =  - ^ -  cos  9  (33) 

al<t)  *  b j_("t) 

distortional  radial  displacement  =  U  (0,  t)  -  aQ(t) - ^ -  cos  6  (34) 

The  average  radial  displacement  is  evaluated  by  numerically  integrating 
the  data  in  figure  23,  to  obtain  the  U  (0,  t)  shown  in  figure  24.  Then, 
aQ(t)  is  evaluated  by  numerically  integrating  the  following  equation 

ao(t)  =¥  [  u  (e,t)de  (35) 

-'0 

and  is  shown  in  figure  25. 

The  rigid-body  radial  translation  is  computed  from  equation  (33)  by 
using  the  values  of  (aj  -  b  j)/2  presented  in  figure  10  of  Reference  6. 
Finally,  the  radial  distortional  displacements  are  calculated  from  equa¬ 
tion  (34)  and  these  results  are  presented  in  figures  26  and  27.  Comparing 
figures  25  and  26,  it  is  noted  that  the  average  and  the  maximum  distor¬ 
tional  radial  displacements,  which  occur  on  0  =  0,  have  about  the  same 
magnitude. 
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Figure  25.  Average  Cavity  Displacement 

(Derived  from  Reference  6  Figure  7  data) 


9'0 
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APPENDIX 

RECIPROCITY  RELATIONS  IN  ELASTIC  AND 
VISCOELASTIC  WAVE  PROPAGATION 
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I.  SUMMARY 


By  a  slight  extension  of  Graffi's  proof  of  the  reciprocity  relations 
in  dynamics,  it  is  pointed  out  that  dynamic  reciprocity  applies  as  well 
to  viscoelastic  media  and  also  to  traveling  loads. 

II.  INTRODUCTION 

The  most  general  reciprocal  theorem  in  dynamics  was  asserted  by 
Horace  Lamb  (Reference  16)  to  be  derivable  from  a  remarkable  formula 
established  by  Lagrange  in  the  "Mecanique  Analytique"  (1809)  by  way  of 
a  prelude  to  his  theory  of  the  variation  of  arbitrary  constants.  Lamb 
showed  how  the  reciprocal  theorems  of  von  Helmholtz  in  the  theory  of 
least  action  in  acoustics  and  optics  and  of  Lord  Rayleigh  in  acoustics, 
can  be  derived  from  Lagrange's  formula.  Rayleigh  (Reference  21)  extends 
the  reciprocal  theorem  to  include  the  action  of  dissipative  forces  and 
Lamb  shows  the  complete  reciprocity  relationship  in  a  moving  fluid  with 
reversed  flow  conditions.  These  theorems  are  best  stated  and  proved  in 
terms  of  generalized  coordinates. 

In  the  theory  of  static  elasticity,  a  well  known  reciprocal  theorem 
is  associated  with  the  names  of  Maxwell,  Betti,  and  Rayleigh,  A  generali¬ 
zation  of  this  theorem  to  dynamic  problems  in  elasticity  is  given  by  Graffi 
(References  12  and  13),  and  certain  applications  of  Graffi's  results  to  the 
problem  of  elastic  wave  propagation,  in  particular,  of  the  reciprocity 
between  certain  results  of  Pekeris'  problem  (Reference  19)  of  a  suddenly 
applied  vertical  load  and  Chao's  problem  (Reference  7)  of  a  suddenly 
applied  horizontal  load,  was  pointed  out  by  Di  Maggio  and  Bleich  (Refer¬ 
ence  9).  However,  there  seems  to  have  been  no  mention  in  the  literature 
about  reciprocity  relations  that  hold  for  traveling  loads.  It  is  the  purpose 
of  this  Appendix  to  state  such  a  relationship,  to  extend  it  to  viscoelastic 
media,  and  to  demonstrate  some  interesting  applications. 
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III.  DERIVATION  OF  EQUATIONS 


Consider  a  linear  isotropic  viscoelastic  medium,  which  is  described 
by  the  stress-strain  relationship  (in  the  usual  tensor  notations) 


3n  ^nJVk 


6..  + 
U 


(1) 


where  cr..  is  the  stress  tensor,  u.  is  the  displacement  vector,  a  ,  b  ,  c 
ij  1  n  n  n 

are  constants  if  the  medium  is  homogeneous,  and  are  functions  of  space 
coordinates  if  the  medium  is  nonhomogeneous.  If  N  =  0,  the  material 
obeys  Hooke's  law  and  the  above  relationship  is  written 


cr. .  =  Xu.  ,6..+u(u.  .  +  u.  .) 

ij  Ac,  k  ij  i,  j  j,  r 


(2) 


where  X,  p  are  the  Lame  constants  which  may  depend  on  space  coordinates. 
If  X.  is  the  body  force  per  unit  mass,  the  equation  of  motion  is 


P 


pX.  +  cr .  . 
i  ij- 


(3) 


If  the  material  occupies  a  finite  volume  V  with  a  boundary  surface 
S  =  S  f  S  ,  surface  traction  being  specified  over  S  and  displacement 

(T  U  O’ 

being  specified  over  S  ,  the  boundary  conditions  that  must  be  satisfied  are 


v 


..n.  =  T. 

on  S 

ij  J  1 

cr 

u.  =  w. 

on  S 

(4) 

i  l 

u 

where  n.  is  the  direction  cosine  of  the  normal  v  of  S.  Let  us  consider 
J  v 

problems  in  which  X ,  (x  ^ ,  x^,  x^;  t) ,  T .  (x  j ,  x^,  x^Jt),  and  w^  (x  ^ ,  x^ ,  x^;  t) 
are  given  functions  of  space  and  time,  v/hich  starts  its  action  at  t  =“0, 
under  the  initial  conditions 


du. 

i 


d  u. 


when  t  —  0 


(5) 
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The  response  m  (x^,  x^,  x^;  t)  is  sought.  To  this  end,  we  apply  Laplace 
transformation  with  respect  to  the  time  t  to  every  dependent  variable 
under  the  assumption  of  suitable  continuity  conditions  so  that  the  trans¬ 
forms  exist.  Let  the  Laplace  transform  of  u  be  written  as  u 


'  oo 

u  =  2^(t)  |  sj=j  e  S  u(t)  dt 

Jo 


(6) 


We  have,  on  account  of  the  initial  conditions  named  above 


where 


-ij 


=  \(s)u,  .  6..  +  u(s)(u.  .  +  u.  .) 

— k,  k  ij  -i,  J  -J,  i 


s  pu.  =  pX.  +  cr. .  . 
— l  —i  -ij ,  j 


cr .  .n .  =  T . 

-ij  J  —i 


u .  =  w . 

—l  — ’i 


in  V 


on  S 


on  S 


u 


7) 


Ms)  =  £  bnsn/^  ansn  ,  P(s)  =  £  cnsn/£ansr 


(8) 


When  the  mathematical  problem  is  solved  for  u.  (x^,  x^,  x^;  s),  the 
inverse  transform  must  be  performed  in  order  to  obtain  the  desired 
solution  u.  (x  j ,  x^,  x^;  t). 


Now  consider  two  problems  in  which  the  applied  body  force  and  the 
surface  tractions  and  displacements  are  differently  specified.  Let  the 
variables  involved  in  these  two  problems  be  distinguished  by  superscripts. 


Then 

.. 

In  V  : 

2  (1)  v(l),  (1) 

s  pu.  =  pX.  +  cr..  . 

-i  —i  -iJ  -  J 

(9a) 

s2pUi( 

2>.  Px.(2>  +  ,r(2>. 

— 1  -1J»J 

(10a) 

On  S  : 
cr 

-IJ  J  -1 

(9b) 

cr^l 
-ij  J 

-  t.(2) 

— i 

(10b) 

On  S  : 
u 

U*1^  w'1* 

—1  — 1 

(9c) 

u.<2>  = 
—1 

(2) 

=  w. 

— i 

(10c) 
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Multiplying  equation  (9a)  by  u.^  and  (10a)  by  u.^,  subtracting,  and 
integrating  over  the  volume  V,  we  obtain 


V 


pX(.  1}u(2)dv  + 
— 1  —l 


<r(.1>.u(2,dv=  pX^u^dvH- 

y-iJ-J'1  JV  _1  'V 


L^.U^ 

J“ 1 


dv  (11) 


Now 


V 


a^.u^dv 

r  •  J-1 


>v 


>,ui1U..)..»!2)tp1!tllu!1>l.««2>ldv 

^k.kij/j-i  \  — l,  j  -J,  1/  J-l  ' 


^U-.u^n-ds-  f  6..u<2>dv 


+  r^u!1!u<2)».d, 

Je  -l.j-l  J  J  c  — J  >  1—1  J 


f  (1)  (2) 
Jv  -1>  J_1>  J 


!dv  -  |  pu^lu^ldv 
y  J  »  1  1  j-  J 


x  (1)  (2)  .  f  .  (1)  (2), 

Xu,  '  u:  n.ds  -  j  Xu,  ,  u:  :dv 

S  “k.k-1  1  Jy  ~ k.k-j.j 


+  f  pu! 1  !ui2^n.ds  +  fpuW.^n.ds 

J  s  1  * J— 1  }  J 

f  pu^lu^dv-  f  pu^lu^dv 

A  similar  expression  is  obtained  for  the  integral 


T^.U^dv 

y-'J'J-1 


(12) 


(13) 


When  these  expressions  are  substituted  into  equation  (11),  we  see  that 
a  number  of  volume  integrals  cancel  each  other.  The  surface  integrals 
contributed  by  the  integrals  ( 12)  and  (13)  to  (11)  are 

f  Xu<1>u<2)n  ds+  fpu*1* u<2)n.ds+  fp uW.^n.ds 
js-k,  k-i  i  Js-i.j-i  J  -J.i-i  J 

=  f  Xu[2U!1)n.ds+  [pu(2?ui1)n.ds+  f  pu<2>u< 1  >n.ds  (14) 
Jc  -k,  k-1  1  J  s  — 1,  J— 1  J  Jq  -J-1-1  j 
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or 


f  (1)  (2)  ,  f  (2)  (1)  , 

tr'.  u:  n.ds  =  I  cr'.  u:  n.ds 

lg  ij  1  j  Jg  Ij  1  J 


(14a) 


If  we  now  recall  the  definitions  of  S  and  S  ,  and  the  boundary  conditions 

cr  u  ’ 

(9b),  (9c),  (10b),  (10c),  and  substituting  all  these  into  equation  (11),  we 
obtain  finally 


pxW.2: 

— l  —l 


dv  + 


(  T^u^ds* 

Jg  1  1 


cr 


Is 


(1)  (2) 

cr: .  w: 

-ij  —l 


u 


n.ds 

J 


f  Pxi2)u*1)dv  +  / 

'  T<2>u(1>ds+  / 

'  rWA.d. 

(15) 

Jv  -1  -1  J 

S  1  1  J 

s  “1J  _1  J 

a  u 


This  is  a  general  reciprocal  relation  in  the  Laplace  transformation  form. 
It  is  of  the  same  form  as  Betti's  reciprocal  relation  in  elastostatics. 

Since  the  inverse  transform  of  the  product  of  two  functions  is  the 
convolution  of  the  inverses,  we  obtain 


u: 


p(x)X^1)(x,  t 


t)u!2^(x,  T)dTdv  + 


Is  J?‘V“ 


T)u!2^(x,  T)dTds 


+  f  [  <t\  } \x,  t  -  t)w!^(x,  T)n.dTds 

JS  Jo  1J  1  J 


f  p(x)x!2^(x,  t  -  t)u!  1  ^(x,  T)drdv  +  f  f  t!2\x,  t  -  t)u!  ^(x,  T)dTds 
JyJo  1  1  Js  Jo  1  1 


+  I  f  \x,  t  -  t)w!  ^(x,  T)n.drds 

Js  Jo  1J  1  J 

u 


(16) 


This  is  the  general  reciprocal  theorem  for  elasto-kinetics.  Whether  the 
material  is  viscoelastic  or  purely  elastic  makes  no  difference  in  the  final 
result.  Note  that  this  result  holds  for  variable  density  p(Xj,  x^,  x^)  and 
nonhomogeneous  material  properties. 
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IV.  GENERALIZATION  TO  INFINITE  REGION 

A  generalization  of  the  above  result  to  an  infinite  or  semi- infinite 
region  is  possible.  Since  a  finite  wave  speed  exists  there  always  exists 
a  boundary  surface  at  any  finite  t  »  0  which  is  yet  uninfluenced  by  the 
loading  initiated  at  t  =  0.  Let  S  be  such  a  surface.  Then  w.  =  0  on  S 
and  the  remainder  of  the  equation  holds  without  question. 


V.  SPECIFICATION  AND  APPLICATIONS 
Space-time  Separable  Body  Forces,  Surface  Tractions,  and  Displacements 
If 


x!u 

l 

=  1 )( 

»g(t) 

x<2> 

1 

=  Z[Z)Mg(t) 

^d) 

i 

=  p!1J( 

;x)g(t) 

^(2) 

i 

=  p(.2)(x)g(t) 

1 

=  w!1* 

l 

(x)g(t) 

(2) 

w. 

l 

bC 

r\J 

V* 

II 

then  equation  (15)  can  be  written,  on  cancelling  If  {_g(t)}  from  every  term,  as 


I  d'PA^dv-,  I  PU,u!2,ds  +  f  w'2V‘»„. 


JV 


L— il  —1 


1  —1 


JS 


1  -lj  j 


ds 


[  p=;!2,u(1|dv  t  (  p!2)u!1,d»  +  f  wl'h^’n.ds 

lv  *  Js  1  Js  *  > 

cr  u 


The  inverse  transformation  gives 


.--(1).  (2), 


P(Zj:  'u:  (x,t)dv+  /  P(  1)u(2)(x,  t)ds  +  I  W{.c,<r\n(x,  t)n.ds  =  lV 

JV  1  1  Js  1  1  Js  1  1J  J 

<7  U 


,(2)  ( 1), 


"\ 


or 


^  pX(1)(x,t)u(2)(x,t)dv  + 
V  1  1 


f  tIVou^I 
.  s  1  1 


x,  t)ds  +  I  w 
S 


(Z)(x,t)a!1)(x,t)n.ds  =  \2 
i  ij  j  lv 
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where  "lN2  indicates  the  same  expression  as  on  the  left-hand  side  except 


that  the  superscripts  (1)  and  (2)  are  interchanged.  Graffi's  well-known 

formula  results  if  wl^,  w|2^  =  0  on  S  . 

11  u 

Forces  Applied  at  Different  Times 

Graffi  also  gives  the  following  case.  If 


x(in  =  z;(11)(x)g(t  -  tl) 


x[2)  =  ^2)(x)g(t  -  T2) 


T(.U  -  p!U(x)g(t  -  T.) 


T(.2)  =  pS2)(x)g(t  -  TJ 


w ; 

l 


m 


(2) 

w:  '  =o 


then  equation  (  1  5)  becomes,  on  cancelling  ^g(t)  |  s^  , 


'  > — ■(  1 )  ST1  (2)  ,  .  f  _(1)  'sTl  (2),  \2 

p -j.  e  u.  dv  t  P.  e  u.  ds  =  1 

v  1  -1  JSa  1 


The  inverse  transform  gives 


(  P,l£(i1)(x)u(i2)(x1t  -  TjJdv  +  J  t  -  Tj)ds 

^  ^cr 


P^’2)(x)u!  ^(x,  t  -  T2)dv  +  J  P^W^fx.t-T^ds 


Concentrated  Forces 

If  the  loading  consists  of  concentrated  loads  f!^  and  fj2^  acting  at 
points  Pj,  p^  respectively,  we  may  consider  or  Pj  as  suitable  delta 

function  and  obtain  at  once 


f(i1)(P1)uS2)(P1-t  -  Ti>  =  f(i2)(Pz)u!1)(p2.,t  -  T2) 
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or,  if  T  j  T2 


fi1)(pi)u|2)(pr  fc)  =  fi2)(p2)u[1>(p2' 


This  is  the  extension  of  the  conventional  elastostatic  Betti-Rayleigh 
reciprocal  relation  to  kinetics. 

Impulsive  and  Traveling  Concentrated  Forces 

Let  an  impulsive  concentrated  force  act  at  a  point  p^ 

px!11  =  fl1)6(p1)6{t) 

and  a  concentrated  force  fj2^  be  applied  at  the  origin  at  t  =  0  and  thereafter 
moved  along  the  axis  at  uniform  speed  U 


pX<2)  =  f[Z)6(t  -  -[ij6(x2)6(x3) 


No  other  surface  loading  or  displacement  is  imposed.  Then  equation  (16) 
give  s 

W..*! 


(1) 


6(x2)6(x3)dx1dx2dx3  6  It  -  -g-  jul  '(Xj ,  x2>  x3>  t  -  T)dT 


and  therefore 


fO)  ,U) 


fl  u-  '(Pl.t)  =  f; 


,(2)  :  U 


J  1 


u.|x,  ,  0,  0,  t  -  |  dx 


u  r  i 


is  known,  then  u!  (p^,  t)  can  be  found  from 


If  u}1)  [xj,  0,  0,  t  -  (Xj/U)]  i 
the  above  equation. 

Suddenly  Started  Line  Load  Over  an  Elastic  Half-Space 


The  steady-state  solution  of  the  problem  of  a  line  load  traveling  at 
constant  speed  over  an  elastic  half-space  (a  two-dimensional  problem) 
was  first  solved  by  Sneddon  (Reference  20)  and  later  by  Cole  and  Huth 
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(Reference  8)  by  means  of  Fourier  transformation  and  by  analytic  functions 
of  a  complex  variable.  A  singularity  was  found  if  the  load  travels  at  the 
Rayleigh  wave  speed  of  the  medium.  To  examine  this  singularity  Ang 
(Reference  1)  considered  the  problem  of  suddenly  started  line  load.  Now 
according  to  the  reciprocal  theorem  Ang's  problem  can  be  solved  by  one 
integration  of  the  solution  of  Lamb's  problem:  the  impulsive  loading  at 
one  point  (not  traveling)  inside  a  two-dimensional  half-space.  Only  the 
surface  displacement  due  to  the  point  loading  needs  to  be  known. 


77777777777  //? 77)/ 7777 

\  fj(2W)  f;(1)S(X-Ut) 

LAMB'S  PROBLEM  ANG'S  PROBLEM 

Figure  1.  Illustration  of  the  Loading  in  Lamb's 
and  Ang's  Problems 

Suddenly  Started  Point  Load  Over  an  Elastic  Half-Space 

The  corresponding  three-dimensional  problem  of  a  concentrated  load 
suddenly  applied  at  a  point  on  the  surface  of  an  elastic  half- space  and 
thereafter  moved  at  a  constant  velocity  U  in  the  x-direction  has  not  yet 
been  solved.  However,  from  the  reciprocal  theorem  named  about,  it 
becomes  apparent  that  the  solution  to  this  problem  can  be  obtained  by  an 
integration  of  Pekeris'  solution  to  a  suddenly  applied  vertical  load  and 
Chao's  solution  to  a  suddenly  applied  horizontal  load  at  a  given  point  in 
an  elastic  half- space. 
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